// --*- C++ -*------x-----------------------------------------------------------
//
//
// Author:          Enrico Negri
//
// Project name:    BiocompUP Align
//
// Date:            03/2008
//
// Description:     A simple program to visualize the score of every position.
//
// -----------------x-----------------------------------------------------------

#include <Alignment.h>
#include <AlignmentBase.h>
#include <GetArg.h>
#include <ScoringP2P.h>
#include <ScoringP2S.h>
#include <ScoringS2S.h>
#include <SequenceData.h>
#include <SubMatrix.h>
#include <iostream>

#include <HenikoffProfile.h>
#include <PSICProfile.h>
#include <SeqDivergenceProfile.h>

#include <LogAverage.h>
#include <CrossProduct.h>
#include <DotPFreq.h>
#include <DotPOdds.h>
#include <Pearson.h>
#include <JensenShannon.h>
#include <EDistance.h>
#include <AtchleyDistance.h>
#include <AtchleyCorrelation.h>


using namespace Biopool;


/// Show command line options and help text.
void
sShowHelp()
{
	cout << "\nALIGNMENT SCORE"
	     << "\nA simple program to visualize the score of every position.\n"
	     << "\nOptions:"
	     << "\n * [--in <name>]   \t Name of input FASTA file"
	     << "\n   [--pro1 <name>] \t Name of target profile (BLAST M6 format) file"
	     << "\n   [--pro2 <name>] \t Name of template profile (BLAST M6 format) file"
	     << "\n   [--out <name>]  \t Name of output file (default = to screen)"
	     << "\n   [-m <name>]     \t Name of substitution matrix file (default = blosum62.dat)"
	     << "\n   [-d <double>]   \t Min. master vs all seq. id. filter threshold (default = 0.00)"
	     << "\n   [-D <double>]   \t Min. all vs all seq. id. filter threshold (default = 0.00)"
	     << "\n   [-u <double>]   \t Max. master vs all seq. id. filter threshold (default = 1.00)"
	     << "\n   [-U <double>]   \t Max. all vs all seq. id. filter threshold (default = 1.00)"
	     << "\n   [--fasta]       \t Use FASTA format to load profiles"
	     << "\n   [--verbose]     \t Verbose mode"
	     << "\n" << endl;
}


int
main(int argc, char **argv)
{
	// Default parameters

	string seq1, seq2, seqA1, seqA2, seq1Name, seq2Name;
	string inputFileName, proFile1Name, proFile2Name, outputFileName;
	string matrixFileName = "blosum62.dat";
	double downs, downa, ups, upa;


	// --------------------------------------------------
	// 0. Treat options
	// --------------------------------------------------

	if (getArg("h", argc, argv))
	{
		sShowHelp();
		return 1;
	}

	getArg("-in", inputFileName, argc, argv, "!");
	getArg("-pro1", proFile1Name, argc, argv, "!");
	getArg("-pro2", proFile2Name, argc, argv, "!");
	getArg("-out", outputFileName, argc, argv, "!");
	getArg("m", matrixFileName, argc, argv, matrixFileName);
	getArg("d", downs, argc, argv, 999.9);
	getArg("D", downa, argc, argv, 999.9);
	getArg("u", ups, argc, argv, 999.9);
	getArg("U", upa, argc, argv, 999.9);
	bool fasta = getArg("-fasta", argc, argv);
	bool verbose = getArg("-verbose", argc, argv);


	// --------------------------------------------------
	// 1. Load data
	// --------------------------------------------------

	string path = getenv("VICTOR_ROOT");
	if (path.length() < 3)
		cout << "Warning: environment variable VICTOR_ROOT is not set." << endl;
	string examplesPath = path + "examples/";
	string dataPath = path + "data/";

	if (inputFileName != "!")
	{
		inputFileName = examplesPath + inputFileName;
		ifstream inputFile(inputFileName.c_str());
		if (!inputFile)
			ERROR("Error opening input FASTA file.", exception);
		Alignment ali;
		ali.loadFasta(inputFile);
		if (ali.size() < 1)
			ERROR("Input FASTA file must contain two sequences.", exception);
		seq1 = Alignment::getPureSequence(ali.getTarget());
		seq2 = Alignment::getPureSequence(ali.getTemplate());
		seqA1 = ali.getTarget();
		seqA2 = ali.getTemplate();
		seq1Name = ali.getTargetName();
		seq2Name = ali.getTemplateName();
	}
	else
		ERROR("aliscore needs input FASTA file.", exception);

	matrixFileName = dataPath + matrixFileName;
	ifstream matrixFile(matrixFileName.c_str());
	if (!matrixFile)
		ERROR("Error opening substitution matrix file.", exception);


	// --------------------------------------------------
	// 2. Output test data
	// --------------------------------------------------

	if (verbose)
	{
		fillLine(cout);
		cout << "\n" << seqA1 << "\n\n" << seqA2
		     << "\n\nRead substitution matrix from file: " << matrixFileName
		     << "\n" << endl;
	}
	fillLine(cout);


	// --------------------------------------------------
	// 3. Select alignment mode
	// --------------------------------------------------

	SubMatrix sub(matrixFile);

	AlignmentData *ad;
	ScoringFunction *sf;
	ScoringScheme *ss;

	if (proFile1Name != "!")
	{
		// Construct first profile

		proFile1Name = examplesPath + proFile1Name;
		ifstream proFile1(proFile1Name.c_str());
		if (!proFile1)
			ERROR("Error opening target profile (BLAST M6 format) file.", exception);

		Alignment ali1;
		if (fasta)
			ali1.loadFasta(proFile1);
		else
			ali1.loadBlastMode6(proFile1);

		if (downs < 999.9)
			ali1.RemoveLowerSimple(downs);
		else
			if (downa < 999.9)
				ali1.RemoveLowerAll(downa);

		if (ups < 999.9)
			ali1.RemoveUpperSimple(ups);
		else
			if (upa < 999.9)
				ali1.RemoveUpperAll(upa);

		Profile *pro1 = new Profile();
		pro1->setProfile(ali1);

		if (proFile2Name != "!")
		{
			// --------------------------------------------------
			// 3.1. Profile-to-Profile case
			// --------------------------------------------------

			// Construct second profile

			proFile2Name = examplesPath + proFile2Name;
			ifstream proFile2(proFile2Name.c_str());
			if (!proFile2)
				ERROR("Error opening template profile (BLAST M6 format) file.", exception);

			Alignment ali2;
			if (fasta)
				ali2.loadFasta(proFile2);
			else
				ali2.loadBlastMode6(proFile2);

			if (downs < 999.9)
				ali2.RemoveLowerSimple(downs);
			else
				if (downa < 999.9)
					ali2.RemoveLowerAll(downa);

			if (ups < 999.9)
				ali2.RemoveUpperSimple(ups);
			else
				if (upa < 999.9)
					ali2.RemoveUpperAll(upa);

			Profile *pro2 = new Profile();
			pro2->setProfile(ali2);

			ad = new SequenceData(2, seq1, seq2, seq1Name, seq2Name);
			sf = new DotPFreq(pro1, pro2);
			ss = new ScoringP2P(&sub, ad, 0, pro1, pro2, sf, 1.00);
		}
		else
		{
			// --------------------------------------------------
			// 3.2. Profile-to-Sequence case
			// --------------------------------------------------

			ad = new SequenceData(2, seq1, seq2, seq1Name, seq2Name);
			ss = new ScoringP2S(&sub, ad, 0, pro1, 1.00);
		}
	}
	else
	{
		// --------------------------------------------------
		// 3.3. Sequence-to-Sequence case
		// --------------------------------------------------

			ad = new SequenceData(2, seq1, seq2, seq1Name, seq2Name);
			ss = new ScoringS2S(&sub, ad, 0, 1.00);
	}


	// --------------------------------------------------
	// 4. Calculate scores
	// --------------------------------------------------

	vector<double> score;

	int i = 1;
	int j = 1;

	for (unsigned int index = 0; index < seqA1.size(); index++)
		if ((seqA1[index] != '-') && (seqA2[index] != '-'))
		{
			score.push_back(ss->scoring(i, j));
			i++;
			j++;
		}
		else
		{
			score.push_back(0.00);
			if (seqA1[index] != '-')
				i++;
			if (seqA2[index] != '-')
				j++;
		}


	// --------------------------------------------------
	// 5. Output scores
	// --------------------------------------------------

	if (outputFileName != "!")
	{
		outputFileName = examplesPath + outputFileName;
		ofstream outputFile(outputFileName.c_str());
		if (!outputFile)
			ERROR("Error opening output file.", exception);
		cout << "\nSaving output to file: " << outputFileName << endl;

		outputFile << " Target   Template               Score\n" << endl;
		for (unsigned int index = 0; index < seqA1.size(); index++)
			if ((seqA1[index] != '-') && (seqA2[index] != '-'))
				outputFile << setw(7)  << seqA1[index]
				           << setw(11) << seqA2[index]
				           << setw(20) << score[index] << endl;
			else
				outputFile << setw(7)  << seqA1[index]
				           << setw(11) << seqA2[index]
				           << setw(20) << "-" << endl;
	}
	else
	{
		cout << "\n Target   Template               Score\n" << endl;
		for (unsigned int index = 0; index < seqA1.size(); index++)
			if ((seqA1[index] != '-') && (seqA2[index] != '-'))
				cout << setw(7)  << seqA1[index]
				     << setw(11) << seqA2[index]
				     << setw(20) << score[index] << endl;
			else
				cout << setw(7)  << seqA1[index]
				     << setw(11) << seqA2[index]
				     << setw(20) << "-" << endl;
	}
	cout << endl;
	fillLine(cout);


	return 0;
}
